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Abstract: Recent analyses of CMB data combined with the measurement of BAO and Hq 
show that dark radiation- parametrized by the apparent number of additional neutrinos 
AN e fj contributing to the cosmic expansion- is bounded from above by about AN e jf < 1.6 
at 95% CL. We consider the mixed axion/neutralino cold dark matter scenario which arises 
in i?-parity conserving supersymmetric (SUSY) models wherein the strong CP problem is 
solved by hadronic axions with a concommitant axion(a)/saxion(s)/axino(a) supermulti- 
plet. Our new results include improved calculations of thermal axion and saxion production 
and include effects of saxion decay to axinos and axions. We show that the above bound on 
AN e ff is easily satisfied if saxions are mainly thermally produced and rriLSP < m-a ^ "is- 
However, if the dominant mechanism of saxion production is through coherent oscillations, 
the CMB data provides a strong bound on saxion production followed by saxion decays to 
axions. Furthermore we show that scenarios with mixed neutralino/axion dark matter are 
highly constrained by combined CMB, BBN and Xe-100 constraints. In particular, super- 
symmetric models with a standard overabundance of neutralino dark matter are excluded 
for all values of the Peccei-Quinn breaking scale. Next generation WIMP direct detection 
experiments may be able to discover or exclude mixed axion-neutralino CDM scenarios 
where s —> aa is the dominant saxion decay mode. 
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1. Introduction 



The recent discovery of a Higgs-like boson with mass rrih — 125 GeV at the CERN LHC 
is an outstanding accomplishment jl], |2j, and seemingly provides the last of the matter 
states encompassed within the Standard Model of particle physics. However, if indeed the 
new boson turns out to be spin-0, then it raises a conundrum: how is it that fundamental 
scalar particles can exist at or around the weak scale when their masses are quadratically 
divergent? The well-known solution is to extend the model to include weak scale super- 
symmetry, which reduces quadratic divergences to merely logarithmic while predicting the 
existence of a panoply of new matter states: the so-called superpartners ||. Indeed, the 
new resonance has its mass sitting squarely within the narrow window m/, ~ 115 — 135 
GeV which is predicted by the minimal supersymmetric Standard Model, or MSSM H|. 

One of the virtues of the MSSM is that it includes several candidates for particle dark 
matter. In i?-parity-conserving theories, it is common to conjecture the lightest neutralino 
Z\ as CDM. Its relic abundance follows from thermal freeze-out in the early universe and 
can be consistent with the observed dark matter abundance, although this imposes severe 
constraints on the MSSM parameter space |J. 

A different problem arises from QCD in that the chiral symmetry U(2)l x U(2)r ~ 
U(2)v x U(2)a associated with two flavors of light quarks gives rise to U(2)v — > SU{2)i x 
U(1)b (isospin and baryon number symmetry) plus a U(2)a whose breaking is expected to 
yield four light pions/pseudo-Goldstone bosons. Weinberg conjectured the U(1)a symme- 
try is violated by quantum effects and indeed 'tHooft showed this was so, explaining why 
rrirj 3> m n . A consequence of 't Hooft's solution to the U(1)a problem is that the QCD 
Lagrangian should contain the CP-violating term 

Ccpv 9 e^-F^FA^ (1.1) 

where g s is the QCD coupling, Fa^v is the gluon field strength tensor (F is its dual) and 6 
is an arbitrary parameter. Measurements of the neutron EDM imply 9 < 10 -10 giving rise 
to the strong CP problem: why is CP violation in the strong sector so small ||? After 35 
years, still the most compelling solution is to invoke an additional Peccei-Quinn symmetry 
U{1)pq whose breaking gives rise to the axion field a J?], || [J: in this case, the offending 
jC-cpv term dynamically relaxes to unobservable levels. Phenomenology dictates the PQ 
symmetry is broken at a scale f a ~ 10 9 — 10 16 GeV, where the lower limit comes from 
astrophysical constraints [ffO] while the upper limit arises from theoretical prejudice | |11| |. 
The axion mass is given by 

/10 12 GeV\ 

m a ~ 6 ^eV . (1.2) 



fa 

Although its mass is predicted to be slight, the axion is still an excellent CDM candi- 
date. It can be produced at temperatures around 1 GeV via coherent oscillations with a 
relic abundance given by [fiX 



f \ 7/6 

^ 2 ^0.23/(g t )^( 10 i2 GeV ) (1-3) 
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where the misalignment angle < 0j < ir and f(6i) is the anharmonicity factor. Visinelli 

r / M 7 /6 

and Gondolo parametrize the latter as f{9%) = In ^ 1 _ggy 7r 2 J • The uncertainty in 
£l s a td h? from vacuum mis-alignment is estimated as plus-or-minus a factor of three. 

In moving towards realistic models, it should be fruitful to incorporate simultaneously 
solutions to both the gauge hierarchy problem and the strong CP problem jnj. In super- 
symmetric axion models, labeled here as the Peccei-Quinn-augmented MSSM or PQMSSM, 
the axion necessarily occurs as but one element of an axion supermultiplet 

a='^^ + iV2ea L + m L T a , (1.4) 
v2 

in 4-component spinor notation Q . Here, we also introduce the -R-parity-even scalar saxion 
field s and the -R-parity-odd spin-1/2 axino field a. In gravity-mediation (as assumed in this 
paper) the saxion is expected to gain mass m s ~ w<3/2> where 777.3/2 is the gravitino mass. 
More generally, the axino mass can range from keV to m 3 / 2 P^L [lq , |1T[] . Here we always 
assume m„ ~ ^3/2 [01) so the lightest susper symmetric particle (LSP) is the neutralino. 
The axion, saxion and axino couplings to matter are all suppressed by the PQ breaking 
scale f a - In the PQMSSM with a neutralino LSP (m a ~ 7773/2 ~ TeV), one expects dark 



matter to be comprised of an axion-neutralino admixture [19, ^y, 21]. In spite of their 
suppressed couplings to matter, both the axino and saxion can play surprising roles in 
dark matter production rates in the early universe. 

In a previous work, Choi et al. showed the effects of axino production on the relic 



neutralino abundance [19]. Thermal production of axinos in the early universe followed by 



their decays at temperatures 

r£ = v / fvWp/(vr 2 ^(r D )/90) 1 / 4 (1.5) 

(where Mp = 2.4 x 10 18 GeV is the reduced Planck mass and T a the axino decay width) 
feeds additional neutralinos into the cosmic plasma if Tg < Tf r , where Tf r is the neutralino 
freeze-out temperature. As a result, the neutralino abundance is augmented with respect to 
the standard thermal (MSSM) abundance. However, if axinos are abundantly produce and 
come to dominate the universe at a temperature T e , the abundance of any relics present at 
the time of axino decay can be diluted by an entropy injection factor r = Sj/So ~ T e /Tg. 
The value of r can range from 1 (no dilution) to values as high as 10 3 — 10 4 or even 
higher |2Cj (see also Ref's ||-|||). In Ref. @, this avenue was pushed much further, 



where analytic formulae were presented for both neutralino and axion production in either 
radiation-, matter- or decay-dominated universes. The effect of saxion production and 
entropy injection from late-time saxion decays was also considered. In Ref. |H|], entropy 
injection from coherent oscillation (CO)- produced saxions with decays s —> gg was found 
to strongly dilute all relics at the time of saxion decay. This effect allowed much higher 
values of f a ~ 10 13 — 10 15 GeV to be cosmologically allowed, even for 9i as large as ~ 0.1. 

While the semi-analytic approach presented in ]l9], ^] and is applicable in many 
cases, there also exist numerous cases where a fully numeric solution to the coupled Boltz- 
mann equations is required. Such cases include the possibility of bino-like neutralinos where 
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(cry) is variable with temperature T instead of constant, and where multiple processes of 
neutralino injection or entropy dilution are possible, such as simultaneously accounting for 
axino, axion, neutralino and saxion production and possible decays. In Ref. pi]] , a cou- 
pled Boltzmann calculation was presented which tracked the abundances of neutralinos, 
axions, saxions, axinos and gravitinos. The calculations were restricted to the case where 
s —> gg and gg were dominant. It was found that SUSY models with a standard thermal 
underabundance of neutralinos (case of higgsino-like or wino-like LSPs) could easily have 
the neutralino abundance enhanced, and in addition any remaining underabundance could 
be filled by relic axions. In such models, the value of f a could be pushed up to 10 13 — 10 15 
GeV owing to the combined effects of neutralino enhancement and dilution from entropy 
production, as well as dilution of relic axions. In the same work, it was found to be difficult 
to suppress the neutralino abundance with respect to the standard thermal abundance in 
the MSSM. At high values of f a where entropy dilution from CO-produced saxions was 
large, a high level of dilution was usually accompanied by a violation of BBN bounds on 
late-decaying relics [32, 34], since it was assumed that most of the saxion energy goes 



into visible energy through the s — > gg, gg decays. However this picture can be significantly 
modified once we consider possible s — > aa, da decays as described below. 



The interactions of saxions with axions and axinos is given by [35] 



\ V PQ , 



+ ••• (1.6) 



with £ = Qi v i / v pq an d where qi is the PQ charge of various PQ multiplets, Vi are their 



vevs after PQ symmetry breaking and vpq = yYli v \q\ = fa/V%- I n some simple models, 
£ can be small or even zero, while in others it can be as high as unity [jig ]. In the case 
where £ is non-zero, there is the possibility of additional decay modes of saxions which can 
influence the dark matter production rate. In particular, for £ > 0.05, the decay s — >■ aa or 
s — > ad can become relevant. The latter decay can feed additional LSPs into the plasma, 
while the former decay gives rise to a population of relativistic axions which forms the 
so-called dark radiation. 

Indeed, up until recently, cosmological data seemed to favor the existence of dark 
radiation, not predicted by the Standard Model. The population of weakly interacting 
relativistic degrees of freedom is parametrized by the number of effective neutrinos (N e ff), 
which is ~ 3 in the Standard Model, corresponding to the three neutrino flavors. Previous 
data from WMAP7, the South Pole Telescope (SPT) and the Atacama Cosmology Tele- 
scope (ACT) suggested N e ff ~ 3.5 — 4.5 [36], indicating a source of dark radiation beyond 



the SM. A variety of papers have recently explored this possibility [[37], |38|, [39, 4C , 35, 41 , ^2| . 
More recently, the ACT jE| has released additional data, which reduced the N e fj value 
combined with the measurement of baryon acoustic oscillations (BAO) and the Hubble 
constant to: 

N eff = 3.50 ± 0.42 (WMAP7+ACT+BAO+# )- (1.7) 



On the other hand, recent SPT [44] and WMAP9 [Eq] analyses reported rather higher 



-3- 



values, 



N eff = 3.71 ± 0.35 (WMAP7+SPT+BAO+# ), (1.8) 
N e fj = 3.84 ± 0.40 (WMAP9+eCMB+BAO+i? )- (1.9) 

From the above numbers it is clear the tension between the latest ACT and SPT/WMAP9 
values for N e ft. While the ACT result has only l.ler- level deviation from the standard 
value, N e fj = 3.04, the SPT and WMAP9 results show almost a 2cr-level deviation 1 . 
Due to the tension between the different analyses and the fact that all current results are 
compatible with the SM model value within 2a, it is hard to consider these results as a 
strong evidence for dark radiation. Instead, we choose to use the current results as upper 
bounds on N e ff, which can be predicted by various models. For the purpose of this paper, 
we will invoke a conservative constraint of 

AN ef f = N eff - Nffi < 1.6 . (1.10) 

Higher values of AN e ff are excluded at 95% CL (or more) by any of the current CMB 
analysis discussed above 2 . 

In this paper, we discuss the impact of the CMB constraint on the PQMSSM parameter 
space, assuming m„ ~ m 3/2 ~ TeV, so the neutralino is the LSP. In order to properly 
compute AN e ff and the cold axion and neutralino relic abundances, we numerically solve 
the coupled Boltzmann equations. Here we make several improvements with respect to the 
results in Ref. 

• we include saxion decays to aa and da final states, 

• we update thermal saxion production rates as recently calculated by Graf and Stef- 
fen and 

• we improve the system of Boltzmann equations to properly compute the amount of 
dark radiation. 

While the decay s — >■ aa may lead to a relativistic component of dark matter which is 
strongly constrained by AN e ff < 1.6, it may also diminish the saxion entropy dilution 
effect which is large when s —> gg dominates instead. Although most of our results are 
weakly dependent on the specific MSSM spectrum, for definiteness we will apply our results 
to two benchmark SUSY models inspired by recent LHC results: one of which contains a 
standard thermal overabundance (SOA) of binodike neutralinos while the other contains 



a standard underabundance (SUA) of higgsinodike neutralinos. As shown in Sec. 5^5, if 



1 It is worth pointing out that in the WMAP9 analysis 'eCMB' denotes the extended CMB, which 
uses the old data sets of SPT (2011) and ACT (2011). Also, each N e ff value is obtained from differ- 
ent data sets for BAO and Ho- Hence it is hard to determine the most updated result. Meanwhile, in 
Ref. ^6), an independent analysis was made, which consistently combines the most recent data sets from 
ACT and SPT with WMAP9 data. The results obtained in this case for ACT+WMAP9+BAO+-Ho and 
SPT+WMAP9+BAO+// are consistent with the latest values reported by ACT Q and SPT (44). 

2 From here on we use 'CMB' to encompass any of the current ACT, SPT or WMAP9 data analysis. 
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s — > aa decays dominate, the SOA case is typically excluded for all f a values by combined 
BBN, dark radiation and dark matter constraints. For the SUA case, the Xe-100 direct 



WIMP search bound [47] also becomes relevant. For SUA, some parameter points evade all 
constraints at low f a ~ 10 9 — 10 11 GeV and also at high f a ~ 10 15 — 10 16 GeV, but otherwise 
the bulk of parameter space points are also excluded by the combined constraints. Next 
generation WIMP direct detection experiments may be able to discover or exclude mixed 
axion-neutralino CDM scenarios with a SUA of neutralinos, if s —> aa is the dominant 
saxion decay mode. 

In the next section we define the two SUSY benchmark points used in our analysis. 
Then, in Sec. [3L we discuss the main saxion and axino decay modes relevant for our results. 
Sec. |] introduces some basic formalism and notation for our subsequent discussion and 
the calculation of AN e ff and the BBN constraints. Analytical expressions for the dark 
radiation energy density in the PQMSSM are derived in Sec. [|, while more general results 
using numerical solutions are presented in Sec. |5.3| . Finally, we present our conclusions in 
Sec. ^. In an Appendix, we provide a detailed description of the formulae for the coupled 
Boltzmann equations used to compute our numerical solutions. 

2. Two SUSY benchmarks 

In this Section, we present two SUSY model benchmark points which are useful for il- 
lustrating the effects of dark radiation: one (labeled as SOA) has a standard thermal 
overabundance of neutralinos, while the other (labeleld as SUA) has a standard thermal 
underabundance. 

For the SOA case, we adopt the mSUGRA/CMSSM model with parameters 

(m , m 1/2 , A , tan/3, sign(fi)) = (3500 GeV, 500 GeV, -7000 GeV, 10, +) (2.1) 

We generate the SUSY model spectra with Isajet 7.83 pq| . As shown in Table [j], the SOA 
point has mg = 1.3 TeV and nig ~ 3.6 TeV, so it is beyond current LHC sparticle search 
constraints. It is also consistent with the LHC Higgs discovery since = 125 GeV. The 
lightest neutralino is mainly bino-like with = 224.1 GeV, and the standard neutralino 



thermal abundance from IsaReD |49| is found to be QM SSM h 2 = 6.8, a factor of ~ 60 
above the measured value [^q| . Some relevant parameters, masses and direct detection 
cross sections are listed in Table |l|. Due to its heavy 3rd generation squark masses and 
large [i parameter, this point has very high electroweak finetuning 



The second point listed as SUA comes from radiative natural SUSY [51] with param- 
eters from the 2-parameter non-universal Higgs model 

(m 0! m 1/2 , A , tan/3) = (7025 GeV, 568.3 GeV, -11426.6 GeV, 8.55) (2.2) 

with input parameters (//, tua) = (150, 1000) GeV. With m g = 1.56 TeV and nig ~ 7 TeV, 
it is also safe from LHC searches. It has = 125 GeV and a higgsino-like neutralino with 
mass = 135.4 GeV and standard thermal abundance Q~ ISSM h 2 = 0.01, low by about 
~ 10 from the measured dark matter density. It has very low electroweak finetuning. 
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SOA (mSUGRA) 


SUA (RNS2) 


m 


3500 


7025 




500 


568.3 




-7000 


-11426.6 


tan /3 


10 


8.55 


t* 


2598.1 


150 


m A 


A OQ A O 


i nnn 
1UUU 


m h 


1 

1ZO 


izo.U 


Uln 
UOg 


1312 


1562 




3612 


7021 




669 


1860 




224.1 


135.4 


Q, s i d h 2 
a SI (Z lP ) pb 


6.8 


0.01 


1.6 x 10~ 12 


1.7 x 10~ 8 



Table 1: Masses and parameters in GeV units for two benchmark points computed with Isajet 7.83 
and using m t = 173.2 GeV. 



3. Saxion and Axino decays 

The partial widths for s — >■ gg and s —> gg have been discussed in several papers (see for 
instance Ref. |]52| ). The decay s — > 77 is also possible but is always subdominant and not 
of consequence here. Instead, we focus on the possibility of s — > aa. Using the Lagrangian 
in Eq. (|1.6| ), we find 

£ 2 mf £ 2 vn? s 
64nv 2 PQ = 327/ 2 



in accord with |35[]. Since it will be relevant for our subsequent discussion, we also list the 



s — )• gg decay width: 



2 3 

a s m s 



^ •"" 32^/,7' (3,2) 
In addition, we can consider the saxion-axino-axino interaction term 

C 3 jL- — sd(])d. (3.3) 
V2wpQ 

From this interaction Lagrangian, we obtain 

r s -,aa = |-^l-4-f . (3.4) 



4vr / 2 V m s/ 

It is worth noting that the sad coupling constant is not necessarily the same as the saa 
one. Indeed, the sad coupling can be generated by both SUSY breaking and PQ symmetry 
breaking so that it is highly model-dependent. The sdd coupling can be much smaller than 



that of saa [53|. In this work, for simplicity, we assume a common coupling constant £, 
since it is sufficient to illuminate the effect of saxion decay to axino pairs. This assumption 
allows us to cover both cases where saxion decays to axion and to axino pairs. 
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Figure 1: Saxion branching fractions s — ► gg (red), s gg (green), s — > aa (dark blue) and 
s — > aa (light blue) versus m s for ma = 0.5 TeV and nig = 1.6 TeV. We take £ = 0.01 (dashed) and 
1 (solid). 



To show the relative importance of these decay modes, we show in Fig. [I] the saxion 
branching fractions vs. m s for the case where = 0.5 TeV and nig = 1.6 TeV. We take 
£ = 0.01 (dashed curves) and £ = 1 (solid curves). In the case of £ = 0.01, we see that 
s —> gg dominates the saxion branching fraction for all m s values. The coupled Boltzmann 



calculation of mixed aZ\ dark matter presented in Ref. [21] pertains to this case. If instead 
£ ~ 1, then the situation changes radically, and we see that s — > aa dominates for all values 
of m s . This decay mode, as mentioned earlier, may lead to substantial dark radiation and 
contribute to N e ff. In addition, once m s > 2rria then s — )• aa turns on and rapidly becomes 
comparable to the s — > aa branching fraction. In this case, saxion decay followed by axino 
cascade decays may feed extra neutralinos into the plasma, thus bolstering the neutralino 
abundance. Since T(s — > aa) ~ while F(s — > aa) ~ m s , then as m s increases well past 
2m„, the s — > aa mode more completely dominates the saxion branching fraction. 

Since here we assume m„ > , the axino is unstable and decays to Z{ + Z/j. Other 
decay channels to charginos and gluinos may also be present if they are kinematically 
allowed. The specific expressions for each of these decay modes have been discussed in 
detail in Ref. ]20[| . Although all of these are included in our results, for most of our 
discussion the only relevant decay mode is a — > g + g with the decay width given by: 

a 2 ( m~ \ 3 

r ~ a ^~™ = T&r% m ^ y 1 ~ J ■ (3,5) 

The above decay is the dominant one as long as it is kinematically allowed. 
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4. Boltzmann equations for energy and number densities 



Here we present a very brief description of our procedure to numerically calculate the relic 
abundance of mixed aZ\ CDM in the PQMSSM and the amount of dark radiation (AN e ff). 
A more detailed discussion is left to Appendix 

4.1 Boltzmann equations 

Under the assumptions described in Appendix |A|, the general Boltzmann equations [54| for 
the number density i%i and energy density pi of a particle species % can be written as: 

^1 + 3Hni = _r f mA + [(ni(T)) 2 - nj]{av)i + V BR(j -> %)T jmj -2- (4.1) 
dt Pi Y Pj 

^ + 3^( Pi + Pi) = (n? - n?) (<rv) — - r^n* + ^ ->• QT^nj (4.2) 

where Pi is the pressure density, Tj the decay width, 7, is the relativistic dilation factor, 
(av) is the (temperature dependent) thermally averaged annihilation cross section times 
velocity for the particle species i, rlj is its equilibrium number density and P>R(j — > i) is 
the branching fraction for particle j to decay to particle i. 3 As discussed in Appendix |A], 
the above equation is also valid for coherent oscillating fields once we take BR(j — >• 1) = 
and {av)i = 0. 

It is also convenient to write an equation for the evolution of entropy (S): 

fry 2 1 \ 1/3 1 

S= [-^9*(T)-) R*J2R(i^X)-T iPi 
R 3 

or S = —Yl R ^ -> X ) T i m i n i ( 4 - 3 ) 

i 

where R(i —> X) is the fraction of energy injected in the thermal bath from i decays. 
We supplement the above with Friedmann's equation: 

1 dR r~PT~ \ - , 7T 2 4 

, with px = ' 1 '' 

p 



" = 5 *- = tfsss • w,th «■= V + m*< t » t ■ (4 - 4) 



where Mp is the reduced Planck mass. The set of coupled differential equations, Eqs. (|4,1|), 
Q4.2| ), (4.3) and (4.4) can be solved as a function of time. At late times (T< 1 MeV), when 
unstable particles have decayed (except for the axion) and the neutralino has decoupled, 
the CDM relic abundance is given by: 

_Pz 1 (T)s(T ) p c a °{T) s{T Q ) 

QcDM — — 7(r)" ~p c ~sir) (45) 



3 In this paper, particle species i denotes 1. neutralinos Z\, 2. thermally produced (TP) axinos a, 3. 
and 4. coherently produced (CO)- and TP-produced saxions s(x), 5. and 6. CO- and TP/decay-produced 
axions a, 7. TP gravitinos G and 8. radiation. We allow for axino decay to gg, yZi and ZZi states 
(i = 1 — 4), and saxion decay to gg, gg, aa and da. Additional model-dependent saxion decays e.g. to hh 
are possible and would modify our results. We assume G decay to all particle-sparticle pairs, and include 
3-body gravitino modes as well |5q]. 
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where p^° is the energy density of cold (coherent oscillating) axions, To is today's temper- 
ature and p c is the critical density. Furthermore, since dark radiation is composed of hot 
axions (both thermally and non-thermally produced), we have: 

A7V Pa(T) p a (T) 120 /11\ 4/3 

where we used p v = ffg-T^ and T v = (jy)^ 3 T. 

4.2 Constraints from BBN and Dark Radiation 

A critical constraint on unstable relics comes from maintaining the success of the stan- 
dard picture of Big Bang nucleosynthesis. Constraints from BBN on late decaying neutral 



particles (labeled X) have been calculated recently by several groups [32, 33, p4|. We 



have constructed digitized fits to the constraints given in Ref. [34] and apply these to late 
decaying gravitinos, axinos and saxions. Typically, unstable neutrals with decay tempera- 
ture below 5 MeV (decaying during or after BBN) and/or large abundances will be more 
likely to destroy the predicted light element abundances. We point out that BBN also 
constrains N e ff since it affects the time of neutron freeze-out and consequently the He/H 
ratio. However this constraint is usually weaker than the CMB constraint we assume here 
(AN eff < 1.6). 

4.3 Sample calculation from benchmark SUA: radiative natural SUSY 

As an example calculation, we adopt the benchmark point SUA from Sec. |2| with a higgsino- 
like neutralino and a standard neutralino underabundance of neutralino dark matter. 

Working in the hadronic axion PQMSSM framework, we assume Tr = 10 6 GeV with 
PQ parameters m 5 = m s = = 2 TeV, 0j = 0.01 and f a = 1.5 x 10 14 GeV. We take 
so = f a , where sq is the initial field amplitude for coherent oscillating saxions. We also take 
£ = 1 so that s —> aa is the dominant saxion decay mode. The various energy densities pi are 
shown in Fig. ^ for % = 7 (radiation), Z\ (neutralinos) , a (combined thermally- and decay- 
produced axions), a co (coherent oscillating axions), s (thermally produced saxions), s co 
(coherent oscillating saxions), a (thermally produced axinos) and G (thermally produced 
gravitinos). The energy densities are plotted against the scale factor ratio R/Rq, where 
Ro is the scale factor at T = Tr. We also plot the temperature T (in GeV) of radiation 
(green-dashed curve). 

We see that at all values of R/Rq the universe in this case is radiation-dominated. 
At T > 1 TeV, the TP axions, saxions and axinos all have similar abundances. At 
these temperatures, the saxion coherent abundance is well above these components while 
the gravitino thermal abundance is far below the other components. At low R/Rq, the 
neutralinos are in thermal equilibrium and their energy density lies well above the abun- 
dances of members of the PQ multiplet. As the universe expands and cools, most com- 
ponents are relativistic, and decrease with the same slope as radiation: pi ~ T~ 4 . The 
exception is the CO-produced saxions, which behave as a non-relativistic fluid and fall- 
off as p c s ° ~ T" 3 . At R/Rq ~ 10 3 , the thermally-produced (TP) axinos, saxions and 
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gravitinos become non-relativistic, so now pT p ^ ~ T 3 . For slightly lower temperatures 
with R/Ro ~ 10 5 , neutralinos begin to freeze-out, and their abundance falls steeply. At 
R/Ro ~ 10 9 , axinos begin to decay and bolster the neutralino abundance. More impor- 
tantly, at R/Ro ~ 10 4 , the decay-produced axion abundance begins to swell due to saxion 
decay until by R/Ro ~ 10 8 the decay-produced axion energy density is nearly equal to the 
radiation density. Also, around R/Ro ~ 10 7 with T ~ 1 GeV, CO production of axions 
begins, and by R/Ro ~ 10 8 , with T <C Aqcd, its energy density begins to fall as T -3 . For 
the case illustrated here, the final neutralino abundance turns out to be 0.12 while the final 
axion density turns out to be 0.01, due to the small value of 9i chosen. The dark radiation 
from decay produced axions turns out to give AN e ff = 0.7 and would be allowed by the 
CMB constraint assumed here. Below we discuss how the dark radiation constrains the 
PQMSSM parameters in more general scenarios. 

5. Dark radiation in the PQMSSM 

As seen in Fig. [l]- unless £ <C 1— saxions mainly decay to axions. Since m s 3> m a and m a < 
meV, the (non-thermally produced) axions injected from saxion decay remain relativistic 
until very late and contribute to N e ff. Furthermore, axions can be thermally produced (in 
equilibrium or out of equilibrium) after reheating and will also contribute as an additional 
relativistic species. Below we discuss under which conditions TP or non-thermally pro- 
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duced (decay-produced) axions lead to an excess in N e ff since such scenarios are severely 
constrained by the new CMB results. First we present approximate analytical expressions 
for AN e ff in order to qualitatively discuss its dependence on the PQ parameters. We then 
use the full set of Boltzmann equations to numerically compute A N e f f and obtain our final 
results. 

When computing the abundance of relativistic axions at late times (T<1 MeV), two 
competing effects must be taken into consideration: axion injection from saxion decays 
and entropy injection from both saxion (s —> gg, gg) and axino decays. While the former 
enhances the amount of dark radiation, the latter dilutes it. Thus, in the case of entropy 
dilution, Eq. fl4.6|) becomes: 

lp a (T) lpJT) 120 4/3 , . 

where r is the entropy dilution factor. The above expression must be computed at T ~ 
eV, which are the temperatures to which the CMB is sensitive. 

As mentioned before, axions can be both thermally and non-thermally produced in the 
early universe, so p a receives contributions from thermal production of axions and thermal 
and coherent production of saxions, followed by s —> aa decays. Thus the axion energy 
density after saxion decays is: 

Z p (T) + BRis -> aa) ( ^H)^ ( 



p a (T) = p^(T) + BR(s^aa)\^^\ ^— 1 Ps (T D ) (5.2) 

where Tb is the saxion decay temperature. In the above expression we do not include the 
possibility of entropy dilution of p a , since this is accounted for by the dilution factor r in 



Eq. (J57T|) . Using 

p s (T D ) = m s Y s s(T D ) and g*s(T ~ eV) = 3.9 , (5.3) 



where s(T) = 27r 2 <7*g(T)T 3 /45 is the entropy density and combining Eqs. (5.1) and (|5.2|), 
we obtain: 

1 n TP 1 / m Y 

AN eff ~ + l%m-BR{s -> aa)g, s (T D r 1/3 ^ (5.4) 

r p u r T D 

The dilution of relics due to entropy injection of saxion and axino decays has been 



extensively discussed in the literature (see Refs. [2C, 31 1 and references therein). In the 
general case both saxions and axinos can dominate the universe and inject entropy at 
different times, which can lead to quite involved scenarios. These will be properly addressed 
once we present our numerical results in Sec. [T^. For our analytical results we assume that 
either saxion or axinos inject entropy (but not both), so we can approximate the entropy 
dilution factor by: 



max 



4 m- a Ya 4 v ^ m s Y s 

1, -R(a — > X) — = — Rls — > X)—— — 

3 V 1 Tg 3 v ' T D 



(5.5) 



and we will later assume that either the axino or the saxion term dominates. In the 
above expression and Y a are the axino decay temperature and yield and R(i — > X) is 
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the fraction of pi which goes into visible energy (photons, leptons and jets). The above 
expression ensures that r > 1, so it is valid even when there is no entropy dilution (r = 1). 
Below we review the analytical expressions necessary to compute p^ p ', Y s and r in order 
to obtain AN e ff using Eq. (pTij). 

The thermal production of axions and saxions (in supersymmetric axion models) was 
calculated in Ref. |35[ | and is given by (for relativistic axions): 

(5.6) 



where s = 2ir 2 5*5 (T)T 3 / '45 is the entropy density. Saxions are also thermally produced at 
the same rate J35|, but since m s 3> m a , saxions become non-relativistic and decay in the 
early universe. Hence their energy density before decay is given by: 

**- - ? ~ - * ***** m (i^v) - m 

In order to compute r, it is also necessary to know the thermal production of axinos, which 
has been computed in Ref. 56 ]: 

^.£ aUxW (l)(^-(_*_)^ ( , 8) 



We point out that Eqs. (5.6)- (|5.q) assume out of equilibrium production of axions, 
saxions and axinos, being valid only if Tr is smaller than the decoupling temperature 
{Tdecj- However, if Tr > Td ec , the axion and saxion number densities are given by their 
thermal equilibrium values: 

^ ~ 5.3 x 10-V5(^) 1/3 r , - - 1.2 x l(T 3 m s and - - 1.8 x l^m a (5.9) 

s s s 

From Eqs. (5.7) and (|5.9|), we can estimate the decoupling temperature: 5 

T dec ^1.4x10" GeV (j^^j ' • (5-10) 

Besides being thermally produced, saxions can also be produced through coherent 
oscillations, resulting in the following energy density: 

Y c ° m - ^ ~ 1 9 x lCT 5 CcY m[n[TR,Ts] ( S ° ) * (5 11) 

Y s m s- g -1.9X10 GeV io8QeV { wl2QeY ) 



4 We stress that the thermal production of saxions and axinos implemented here are calculated in different 
ways. The saxion yield from Ref. |5q| is obtained from the Hard Thermal Loop (HTL) approximation 
method, while the axino yield from Ref. |5(| is obtained from finite temperature field theory. According 
to Ref. g|, the HTL yield may be smaller than the finite-temperature calculation by factors as large as 
~ 4 for low reheat temperatures (Tr ~ 10 6 GeV), when the strong coupling constant g s ~ 1 and the 
HTL calculation is no longer valid. Since a finite-temperature calculation for the axion/saxion yield is not 
presently available, we use the result of Ref. | p5[ , even for low Tr values. 

5 Due to the different calculation methods used to compute the non-thermal production of saxion/axions 
and axinos mentioned above, axinos decouple at slightly smaller temperatures than axions and saxions, but 
for simplicity here we take a common decoupling temperature for axions, saxions and axinos. 
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where T s is the temperature at which saxions start to oscillate, given by 3H(T S ) = m s . In 
the above expression, sq is the initial saxion field amplitude and depends on the UV details 
of the model and the inflation dynamics. Natural scales for so are usually taken to be of 
order the PQ breaking scale f a or the reduced Planck mass Mp |5?J . 

Finally, to compute AN e ff we also need the axino and saxion decay temperatures. 
Assuming that these fields decay in a radiation dominated universe 6 , we have: 

Tg = v'PaMp/^VCTlO/gO) 1 / 4 and T D = Vr s M P /(vr 2 ^(T D )/90) 1 / 4 (5.12) 

where r„ and T s are the axino and saxion widths, respectively. 

Using the expressions presented above, we can compute the amount of dark radiation 
in PQMSSM models: 

AN eff ~ ANjfi + 18.02-BR(s ->■ aa)^^)' 1 ^ ^^ + Y ^ (5.13) 

it Id 

where ANj P j = ~Pa P /Pi>, Y^° and Y<f p are the coherent oscillation and thermal yields 
of saxions and 



max 



' 3 v 1 T% 3 v ' T D 



(5.14) 



Before discussing the dark radiation constraints to the PQMSSM parameter space, we point 
out that the first term in Eq. ( |5.13 ) is always subdominant. The maximum contribution 



of thermally produced axions happens when these are produced in equilibrium (Pa P = p a ) 
and there is no entropy dilution (r = 1), which results in the following upper bound: 

^<^00 4/3 ^*lO- (5.15) 



where we used Eq. ( |5.9| ) with 5*5 (T ~ eV) = 3.9. As we can see, the contribution from 
TP axions is well below the current experimental sensitivity and can be safely neglected 
in our subsequent discussion. Hence, in the results below, we only consider the saxion 
contribution to AN e ff. 

1 , m (Y co -4- Y TP ) 

AN eff ^ l8.02-BR(s -> aa)g* 3 {T D )-V 3 Tg • (5-16) 

In most cases, AN e ff is dominated either by Yp° or Yj F so it is interesting to 



separately discuss each of these scenarios. In Sec. |5.1| , we discuss the case where the main 
contribution to AN e ff comes from thermal production of saxions (Y s ~ Y pp ), while in 
Sec. 5.2 we discuss the case where Y s ~ Yg . Then, in Sec. p\l, we present our numerical 



results which do not assume the sudden decay approximation, properly take into account 
the effects of axino and/or saxion dominated universes and cover all possible scenarios, 
including Y^° ~ Y pp . 



6 It is also possible that saxions and/or axinos decay in an axino or saxion dominated universe, which 
modify the values of To and T D . However, in this section we neglect these effects. Once we present our 



numerical results in Sec. 5.3, these effects are automatically taken into account. 
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5.1 AN e ff from Thermal Production 



Here we assume that saxions are mostly thermally produced, so Y^F P 3> Y^° . This happens 
for large Tr and/or small sq. In this case, Eq. ( 5.16 ) becomes: 



1 



AN eff ^ l8.02-BR(s -»• aa)g* s (T D ) 



l/3 /,t s J s 



(5.17) 



We also assume that £ > 0.05, so saxion decays to gluons and gluinos are suppressed (see 
Fig. pi). In this case, entropy injection is dominated by axino decays, so that 



max 



1, 



3 Tl 



(5.18) 



where we have assumed Ria — > X) ~ 1 which is usually a good approximation unless 
m a — Trig ■ 

We will show below that thermal production of saxions usually gives AN e ff < 1, so 
the CMB constraint is easily satisfied in this case. In order to show this, it is sufficient to 
compute an upper bound for AiV e //. From Eqs. ( 5.17 ) and ( 5.18j ) we have: 



AiVef/ < l8.02BR(s -> aa)g* s (T D ) 



l/3 m s Y s TP ,Am~ a Yl p 



/( 



3 T« 



(5.19) 



Note that the equality is satisfied in the case of entropy dilution (r > 1), while it overes- 
timates AN e ff if r = 1. Now, from the expressions for the thermal production of saxions 
and axinos, Eqs. ([Tt|), (tEH) an d (|5.9[) we have: 



yF/yF < Y s /Y- a 



(5.20) 



where is the thermal equilibrium yield for i. Hence: 

AN eff < 9MBR(s -> aa)^ s (rD)- 1/3 — ^ 

ma Id 



9MBR(s -> ao)^5(2b) _1/1 Vs(2B 



-ION -1/4 



"In 




(5.21) 



where we have used Eq. ( 5.12j ) for the axino and saxion decay temperatures. In order to 
simplify even further the above expression, we just need to compute T a /T s . From Fig. [j] 
we see that the s — > gg is always subdominant, while both s — > aa and s — >■ gg can be 
dominant, depending on the value of £. Furthermore, unless m s ~ 2m a , the decay to axinos 
can also be neglected, which we do here for simplicity. Hence we can approximate the total 
saxion width by: 



m" 



327T/2 V ^ 

which gives the following branching ratio for the decay to axions: 

e 



BR(s — > aa) 



(5.22) 



(5.23) 
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As discussed in Sec. g the axino decay width depends on the gaugino spectrum and 
the axino mass. Several axino decay modes have been computed in Ref. pGf | where it has 
been shown that the axino decay width is dominated by a — > gg, if m„ S> m g . Thus, if we 
take the limit m g — > 0, it is easy obtain an upper bound for T„: 



o 



2 



Fa < r- a ^ g g(mg "> 0) = T^f^ ^1 • (5.24) 



16vr 3 / ( 2 

Finally, using the above results in Eq. (|5.21|), we have: 



AAr„,< (0.09 - 0.19) (e+ ^ 2f ^ ^ (5.25, 

where the above range corresponds to g*s(TD), g*s{Tp) = 10 — 100. In Eq.( |5.25| ), the 
equality is satisfied only if Tr > Td ec , r > 1 and m„ 3> m g . Otherwise, AN e ff is below 
the quoted value. Nonetheless the above result illustrates the fact that for models with 
£ ~ 1 and m„ < m s , the CMB constraints are automatically satisfied. It is interesting 
to notice that, although smaller values of £ suppress BR(s —> aa), the amount of dark 
radiation actually increases, due to the increase in the saxion lifetime, which compensates 
the decrease on the branching ratio. We also point out that although the above bound 
is quite general and independent of all other PQ paramaters, it is only valid for thermal 
production of saxions. As we will see in the next section, this result drastically changes 
once we consider coherent production of saxions. 

In order to verify the above results, we numerically compute AN e jf as a function of f a 
using the coupled Boltzmann equations for the PQ fields, but neglecting the contribution 
from CO saxions. We take m s = 1 TeV, Tr = 10 10 GeV, £ = 1 and m„ = 3 and 32 TeV. In 
Fig. ^, the solid blue (red) line shows the full numerical solution for m„ = 3 (32) TeV, while 
the dashed lines show the respective upper bound from Eq. ( 5.25 ) using the appropriate 
values of g*s(Tr)) and (/^(Tg). The dashed gray line shows the amount of entropy dilution 
from axino decays for m„ = 3 TeV. As we can see, AN e ff is always below the upper 
limit from Eq. ( 5.25| ) even when the full numerical solution is considered. Furthermore, for 
most values of f a it is well below the bound since in these regions we have r = 1 and/or 
Tr < Trf ec , where Eq. (|5.25| ) is too conservative. 

The results from Eq. ( |5.25| ) and Fig. || show that in order to violate the CMB constraint 
on dark radiation (AN e fj < 1.6), the axino needs to be at least two orders of magnitude 
heavier than the saxion, if £ = 1. This is hard to achieve on most supersymmetric PQ 
models, where typically mj < m 3 / 2 ~ rn s [jl?]]. For smaller £ values, AN e ff increases but 
it is still below the CMB bound as long as ma < m s . To illustrate this we show in Fig. [| 
the maximum allowed value of rria/m s as a function of £ according to the analytical result 
of Eq. ( |5.25 ). The region below the curve satisfies AN e ff < 1.6, irrespective of the other 
PQ parameter values, as indicated by Eq. ( 5.25| ). On the other hand the region above the 
curve can be either allowed or excluded depending on the choice of PQ parameters. As 
shown in Fig. |], the CMB constraint can be easily satisified for any value of £ if m« < 2m s . 
We also point out that if the current ACT j£J allowed interval for N e ff (AN e ff < 1.3) is 
assumed, it can still be satisfied for any value of £ as long as m„ < 1.4m s . Therefore, we 
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Figure 3: AN e ff as a function of the PQ breaking scale f a for thermal production of saxions only. 
The fixed PQ parameters are m s = 1 TeV, T R = 10 10 GeV and £ = 1. The solid blue (red) line 
corresponds to the numerical solution for ma = 3 (32) TeV. The dashed lines correspond to the 
respective analytical upper bounds given by Eq. (5.25). The entropy dilution factor (r) for ma = 3 
TeV is shown by the dotted gray line. The light blue region has N e ff > 1.6 and is excluded at 95% 
C.L. by the CMB results. 



conclude that the CMB constraint on AN e ff can be easily accomodated in the PQMSSM 
if saxions are mainly thermally produced. In the following we discuss the case in which 
saxion production is dominated by its coherent oscillation component. 

5.2 AN e ff from Coherent Oscillations 

As shown in the previous section, the contribution to AN e ff from thermal production of 
axions and saxions is suppressed unless m„ 3> m s and/or £ ~ 0.05. In this section, we 
assume m„ < m s and £ ~ 1, so AN e ff from the thermal production is negligible and 
the relic density of relativistic axions is dominated by coherent production of saxions and 
their decay. Furthermore, we consider the case for <C Yf to exclude the dilution 
from axino decay. However, even if the axino thermal production is suppressed, axinos can 
still be non-thermally produced through saxion decays, if m s > 2ma. In this case both 
CO saxions and non-thermally produced axinos may inject entropy in the early universe 
at different times, what considerably complicates the picture. Such scenarios will be fully 



addressed once we present our numerical results in Sec. 5^. Here, for simplicity, we will 
assume m s < 2rria, so both the thermal and non-thermal production of axinos can be safely 
neglected. At the end of this section we will briefly discuss what happens if m s > 2mj. 
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Figure 4: The upper bound on mj/m s as a function of £ according to Eq. (5.25). The shaded 
region below the curve satisfies AN e tf < 1 .6, irrespective of the other PQ parameter values. The 
region above the curve may be either excluded or allowed by CMB constraints depending on the 
values of Tr, f a , m s and trig. 



Under the above assumptions Eq. ( |5.16D becomes: 

1 m Y co 

AN eff ~ l8.02-BR(s -> aa)g^s{T D )- X ^ ° s 
r T D 



(5.26) 



and 



max 



'3 1 ' T D 



(5.27) 



Since the only invisible decay mode of the saxion is s — > aa, we can approximate the 
fraction of visible energy injected from saxion decays by: 



R(s->X) = 1- BR(s 



aa) 



e 



a 2 /n 2 



£ 2 + aj/n 2 £ 2 + a 2 /Tr 2 



(5.28) 



where we used the result from Eq. ( 5.23 ) and have once again neglected the decay into 
axinos, since here we assume m s < 2toq. Combining Eqs. ( 5.26 ) and ( 5.27] ) we obtain: 



AiV e// ~ 18.02^5 (T D )- 1/3 < 



4 l ^F 1 



, if r = 1 
, if r > 1 



(5.29) 



Since ff*s(Tb) -1 / 3 > 0.1 and q. 2 /tt 2 ~ 10 -3 , once s ^ aa dominates over the visible mode 
s — > gg, i.e. £ > 0.05, the above result shows that, for coherent production of saxions, 
the case of r > 1 is automatically excluded by the CMB constraint. This means that the 
universe can never have a saxion dominated era and has always been radiation dominated 
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until late times (T<1 MeV). We stress that this conclusion holds as long as saxions are 
mainly produced through coherent oscillations (Y^ p <C Y s ) and m s < 2m a . 

On the other hand, if r = 1 (no saxion dominated era and no entropy injection), the 
constraints on AN e ff give an upper bound on Y^ m s /Ti). In order to rewrite this in 
terms of the PQ parameters, we use Eqs. ( 5.11| ), ( |5.12| ) and fl5.22| ), which gives: 



m»Y„ k 1 xi /a /min[Tf?, TJ \ 

-55- * 2 - 2 x 10 "7FTl7^* 5(rD) (^w) 

10 3 GeV\ 3/2 / s x2 / / a \ 

V (lo^Gev) Vlo^GevJ (5 - 30) 



m. 



where T s ~ 1.3 X lO^m^/lO 3 GeV is the saxion oscillation temperature [p2[ . Thus, from 
Eq. (|5T29D for r = 1: 



k. £ 2 / min[Tr?,T, 

10 3 GeV\ 3/2 / s \V fa 



m s J V10 12 GeV/ V!0 12 GeV 



(5.31) 



where we took <7*s(1d) ~ 100. 

To illustrate the above results and check their validity, we show in Fig. || the numerical 
solution for AN e ff computed using the full set of coupled Boltzmann equations. We assume 
m s = 2 TeV, m- a = 1.5 TeV, m- g = 1.6 TeV, T R = 10 6 GeV, s = f a = 3 x 10 14 GeV, 
6i = 0.01 and vary £. We take the SUA benchmark point. The dotted gray line shows the 
entropy dilution factor, r, and the dashed lines the analytical solutions from Eqs. ( 5.31| ) 



and ( ggg ). As we can see, for the high f a chosen, AN e ff violates the CMB constraint 
for £ > 0.02. Also we see that the numerical solution follows the behavior expected from 



Eq. (5.29) for r > 1 (£ < 0.1). For higher values of £, the s — > gg is extremely suppressed, so 



there is no entropy injection and the solution follows the behavior described by Eq. ( 5.31| ) 



instead. The small difference between the analytical and numerical solutions is expected, 
due to the approximations used in the analytical calculation. We also show the values of 
the neutralino and CO axion relic densities. In this case, since the s — > da and s — ^ gg 
decays are kinematically forbidden, there is no neutralino injection from saxion decays and 
its relic density depends on £ only through the entropy dilution factor. As we can see both 
h 2 an d ^Zi^ 2 SC8 ^ e as V r > as expected. 

The impact of the CMB constraints on the PQ parameter space for £ = 1 is summarized 
in Fig. |6l where we show the excluded region in the / a -?ris plane for different values of Tr 
and sq- The constrained region is computed using the analytical results of Eq. ( |5.31 ). For 



s o = fa and m s ~ 1 TeV, the CMB constraint on dark radiation requires f a < 10 12 — 10 14 
GeV depending on Tr. On the other hand, if the saxion field amplitude so is ~ Mp/100 (as 
suggested in some models |3C] and assumed in Fig. ^6) we see that relativistic axions from 



saxion decays easily violate the CMB constraint, unless m s > 200 GeV, f a < 10 13 GeV 
and Tr < 10 6 GeV. As we can see, the constraint AN e fj < 1.6 imposes an upper bound on 
f a , which strongly depends on the value of sq, since this parameter controls the amplitude 
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Figure 5: Numerical solutions for AN e ff, h and as a function of £ for CO-produced 

saxion case. The dotted gray lines show the entropy dilution factor, r, while the dashed lines show 
the analytical results from by Eqs. (5.31) and (5.29). We consider the parameters m s — 2 TeV, 
m~ a = 1.5 TeV, rrig = 1.6 TeV, T R = 10 6 GeV, s = /, = 3x 10 14 GeV and 0; = 0.01. The MSSM 
point assumed is the SUA listed in Table |l|. 



of saxion coherent oscillations. We point out that this constraint is independent of the 
misalignment angle &i and is not related to the traditional upper bound on f a (< 6^ 2 10 12 
GeV) coming from the overclosure of the universe from CO production of axions. Thus 
the dark radiation constraint provides an additional and independent constraint on f a . We 
also point out that, although in Fig. || we have assumed £ = 1, we expect even stronger 
constraints for £ < 1, since AN e ff increases as £ decreases (as long as £ > 0.05), as shown 
in Fig. H 

So far all the results presented in this section have assumed m s < 2ma- As mentioned 
before, if the s — > da decay is kinematically allowed, saxion decays to axinos can result in 
an even later entropy and neutralino injections (through the axino cascade decay), even 
if thermal production of axinos is suppressed. Since in this case there are two phases of 
entropy injection (at the saxion and axino decays) and the universe can go through saxion 
and axino dominated eras, it becomes difficult to treat it analytically. Therefore, to discuss 
the m s > 2mj case we use the numerical method discussed in the Appendix. We show 
in Fig. |^a the numerical solutions for AN e ff, ft^°h 2 , h 2 and r as functions of m s . 
We take m- g = 1.6 TeV, T R = 10 6 GeV, s = f a = 10 15 GeV, 0< = 0.01 and £ = 1. 
For m s < rria we fall into the scenario discussed above and since we have r > 1 in this 
region, AN e ff is approximately constant, as expected from Eq. ( |5.2£ ). Also AN e fj » 1, 
as anticipated by our previous results for the case r > 1. However, once m s > 2ma, the 
s — > da decay becomes kinematically allowed and drastically reduces AN e ff. This is mainly 
due to two reasons. First, as seen in Fig. [j], BR(s — > aa) decreases around m s > 2rria, thus 
decreasing the injection of axions from saxion decays. Second, the injection of axinos and 
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Figure 6: Dark radiation bounds on the f a -m s plane for so = f a (left) and so = Mp/100 (right). 
The region above the solid red (blue) line has AN eff > 1.6 for T R = 10 10 GeV (T R = 10 6 GeV). 
The curves assume m s < 2rria, only include the contribution from CO saxions and were obtained 
using Eq. (5.31) and £ = 1. 



their subsequent cascade decay to neutralinos significantly injects entropy at later times, 
thus suppressing AN e ff. The sudden enhancement in entropy injection once s — > da opens 
up is shown by the rapid increase in r around m s = 1 TeV. 

To illustrate these effects, we show in Fig. [7|b the cosmological evolution of the energy 
densities of radiation, neutralinos, axions, saxions and axinos as a function of the scale 
factor R. The PQ parameters are the same used in Fig. 0a and m s = 2 TeV. As we can 
see, as saxions start to decay (around R/Ro ~ 10 3 ), the energy density of axions and axinos 
rapidly increases until R/Ro ~ 10 9 , where the universe goes from a saxion dominated to 
an axino dominated era. At much later times (R/Rq ~ 10 13 ) the axino decays and injects 
entropy, significantly diluting the relic density of both relativistic and cold CO axions. As a 
result, AN e ff and Q^°h 2 become highly suppressed, as seen on Fig. [?]a. Therefore, due to 
the late entropy injection from axino decays, the CMB constraint on AN e ff can be easily 
satisfied even for a f a value well above the bounds from Fig. ^. However, the injection 
of neutralinos from non-thermally produced axinos easily surpass the observed DM relic 
abundance, as shown in Fig. [7|. As a consequence, it seems difficult to simultaneously 
satisfy the CMB constraints on AN e fj and Q,£)Mh 2 for large f a . In the next section we 
will generalize this result for more arbitrary choices of the PQ parameters. 

5.3 Numerical Results 

In the previous sections, we have discussed some particular scenarios and derived useful 
analytical approximations for AN e ff. Here we generalize these results using the full numer- 
ical solutions for the coupled Boltzmann equations, simultaneously including all production 
mechanisms (thermal scatterings and coherent oscillations) for the PQ fields. We will also 
discuss the DM content of the viable scenarios. To keep our results as general as possible, 
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Figure 7: Left: AN eff , n^°h 2 , ft^h 2 and r as functions of m s for the CO-produced saxion case. 
We consider the PQ parameters T R = 10 6 GeV, s = f a = 10 15 GeV and 6, = 0.01 and the SUA 
point. Right: energy densities versus the scale factor for the same PQ parameters and m s = 2 TeV. 

we scan over the following PQ parameter values: 

10 9 GeV < f a < 10 16 GeV, 
0.3 TeV < m- a < 20 TeV, 
0.3 TeV < m s < 20 TeV, 

10~ 4 < S / f a < 10 4 , 
10 6 GeV < T R < min(/ a , 10 10 GeV) (5.32) 

while keeping £ = 1 and mg = 3 TeV fixed as well as the MSSM spectrum. For each point 
in parameter space, we use the numerical solutions of the Boltzmann equations to compute 
AN e ff and other quantities of interest. 

In Fig. [8| we once again show AN e ff vs. f a , but now varying all the PQ parameters 
within the parameter space denned in Eq. ( 5.32j ) and using the SUA parameters for MSSM 



spectrum. In order to compare the full results to the analytical approximations of the 
previous sections, we show in different colors points with sq/ f a = 10~ 4 (blue), sq// < 1 
(magenta) and sq/ f a > 1 (red). For the blue points, the low sq/ f a value suppresses the 
coherent production of saxions, so the main contribution to AN e ff comes from thermal 
production, except for very high f a values (> 10 14 GeV). These points are described by 



the results from Sec. 5.1 and the blue dashed line shows the maximum AN e tf allowed by 
Eq. (|5.25| ) (~ 0.8, for rria/m s < 20/0.3 and £ = 1). As we can see, this upper bound is 
respected by the numerical solutions even when all the PQ parameters are varied, except 
for the cases where there is a signficant contribution from coherent production of saxions 
(magenta and red points). 



From Fig. ||, we see that the conclusion from Sec. 5.1 is preserved even in a more general 
scan: thermal production of saxions can easily satisfy the CMB constraint. However, once 
production of saxions through coherent oscillations is included (represented by magenta 
and red points at low f a and by all points at large /„), large values of AN e ff can be 
generated. Nonetheless we find that for sufficiently low Tr and heavy saxions, values of f a 
as large as 10 16 GeV can still be consistent with the CMB constraint even for sq > /„. All 
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Figure 8: AN e ff as a function of the PQ breaking scale f a for the scan over the PQ parameter 
space defined in Eq. ( [5.32 ) for the SUA benchmark. Blue points have So/f a = 10~ 4 and are 
dominated by thermal saxion production for f a < 10 14 GeV. Purple points have so/f a < 1 and 
red points have sq/ f a > 1- The dashed blue line shows the maximum expected value for AN e ff 
from thermal saxion production, as obtained from Eq. ( 5.25| ). The shaded region violates the CMB 
constraint on dark radiation. 



these solutions have m s > 2rria and correpond to the cases where entropy injection from 
axino decays highly suppress AN e ff, as shown by the examples in Fig. f^. 

In Fig. |9], we show the neutralino relic density as a function of f a again for the bench- 
mark point SUA, which has an standard underabundance of neutralino DM due to a 
higgsino-like neutralino. Blue (red) points are allowed (excluded) by BBN constraints and 
have AN e ff < 1.6, while magenta points have AN e ff > 1.6. The green points are both 
allowed by BBN and lie in the la interval for AN e ff from the current WMAP9 results. 
The standard thermal value for Qgji 2 is shown by the dashed gray line and we see that 
f° r fa ^ 10 13 GeV the neutralino relic abundance is enhanced by TP axino decays, s — > aa 
and/or s — >■ gg decays. For larger values of f a , there are several solutions with suppressed 
values of ^^/i 2 when compared to the MSSM value. These points usually have suppressed 
axino and thermal saxion production (due to the large f a value) and s — > gg is forbidden 
(m s < 2rrig). In this case, the injection of neutralinos from axino and saxion decays is 
highly suppressed and easily compensated by the entropy injection from CO-produced sax- 
ions followed by decays to gluons. However, as shown in Fig. [], all these points have too 
large values of AN e ff. This is in agreement with the results of Sec. [5]^, where we showed 
that it is not possible to have entropy dilution (r > 1) from coherent oscillating saxions 
without either violating the CMB constraint on dark radiation or overdosing the universe 
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Figure 9: f2g i h 2 as a function of the PQ breaking scale f a for the scan over the PQ parameter 



space defined in Eq. ( fc.32| ), 

assuming the benchmark point SUA, listed in Table |]. Blue and red 
points have AN e ff < 1.6, while green points have 0.4 < AN e ff < 1.2 and magenta points have 
AN e ff > 0.1.6. Also, blue and green points are allowed by the BBN constraints on decaying 
saxions, axinos and gravitinos, while red points are excluded. The gray dashed line shows the 
standard thermal value fl% p h 2 in the MSSM. The blue-shaded region is excluded by Xe-100 WIMP 
searches at mj? = 135.4 GeV after applying a re-scaled local WIMP abundance. 



(% h 2 > 0.11). The blue-shaded region in the Figure is excluded by applying the recent 



Xe-100 WIMP search bounds J47| to SUA with a re-scaled local abundance of WIMPs. As 
we can see, for the SUA point, the large annihilation cross-section required to suppress 
the neutralino relic is related to a high a SI (Z\p), hence this point is subject to stringent 
bounds from Xe-100. 

Therefore, we conclude that in order to have AN e ff < 1.6, the neutralino abundance 
can only be enhanced with respect to its thermal value. Hence only SUSY models with 
a standard underabundance of relic neutralinos can be consistent with the simultaneous 
constraints on the dark matter density, AN e ff and the constraints from BBN. This result 



is confirmed by Fig. 1C where we again show the neutralino relic density as a function of 
f a , but now for the SOA benchmark, which has a bino-like neutralino with a standard 
overabundance (nf SSM h 2 = 6.8). As we can see, all points consistent with the observed 
CDM abundance are excluded by the AN e jf < 1.6 constraint. 

Finally, we comment on the possibility of the excess in N e ff seen by the WMAP9 and 
SPT analyses being real. In this CclSG, clS shown by the results from Sees. |0| and such 
an excess can only be generated by CO-production of saxions, unless m a 3> m s . Also, as 
shown by the results in Fig. 0, AN e ff ~ 1 can only be obtained if the saxion decay to 
axinos is suppressed (m s < 2rng), otherwise entropy injection from NTP axinos efficiently 
dilute AN e ff to extremely low values. Furthermore, as shown by Fig. ||, AN e ff ~ 1 can 
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Figure 10: Same as in Fig. ^ but for the benchmark point SOA. 



be obtained for a wide range of f a values, depending on the value of sq. This is also seen in 
Fig. |jj where we show as green points the solutions which simultaneously satisfy the BBN 
constraints and lie in the la interval for AN e ff given by the WMAP9 analysis. As we can 
see, AN e jf ~ 1 can be easily obtained as long as f a > 10 10 GeV, where CO-production 
of saxions becomes sufficiently large to generate the excess. We point out that this lower 
bound on f a is directly related to the maximum value of sq (< 10 4 / a ) and the £ value 
assumed in our scan. For higher values of sq and smaller values of £ (< 1), lower f a values 
would be consistent with AN e ff ~ 1. From Fig. |9]we also see that the allowed range of f a is 
more strongly constrained by DM and XelOO bounds than by the AN e ff ~ 1 requirement. 
Therefore we conclude that, if AN e ff ~ 1 is eventually confirmed by Planck data and 
< nia < m s is assumed, saxions were coherently produced in the early universe at 
large rates and m s < 2ms. 



6. Conclusions 

In this paper, we have presented the impact of the CMB constraints on dark radiation 
for the mixed neutralino/axion dark matter scenario. To properly compute the dark mat- 
ter and dark radiation relic abundances we have made use of a simplified set of coupled 
Boltzmann equations, as described in the Appendix. 

We discussed the case of large £ (> 0.05) where saxion decays to aa and da are 
dominant. The case of small £, where saxions dominantly decay to gg or gg was previously 
presented in Ref. [21]. In the present case, s —> aa may contribute to dark radiation 
with a contribution parametrized as AN e ff, the non-standard contribution to the number 
of effective neutrinos. Recent CMB analyses restrict AN e ff < 1.6, providing a strong 
constraint on models producing dark radiation. 



- 24 - 



Our main results may be summarized as follows. 

• At low f a and so (~ 10 9 — 10 11 GeV), saxion production is expected to be dominated 
by thermal production. In this case, the AN e fj < 1.6 bound only weakly constrains 
the PQMSSM parameter space as summarized in Fig. [|. 

• At high f a and so (~ 10 12 — 10 16 GeV), thermal production of saxions, axinos and 
axions are suppressed and coherent production of saxions dominates. In this case, 
the constraints on dark radiation provide an upper bound for f a and sq, as shown in 
Fig. |. 

Once the CMB constraints on AN e ff are combined with the BBN constraints on late 
decaying particles and the constraints on the dark matter relic abundance, we have found 
that: 

• In the case of SUSY models with a standard overabundance of neutralinos (as in our 
SOA benchmark case), it is possible to dilute the relic abundance of neutralinos below 
the observed DM relic abundance through entropy injection from saxion decays to 
gluons. However such solutions always violate the CMB bound on AN e ff. Therefore 
we find that SUSY models with a SOA of neutralinos are still excluded for all choices 
of PQ parameters. 

• In the case of SUSY models with a SUA of neutralinos, low values of f a (~ 10 9 — 10 12 ) 
can easily accommodate all the constraints, since in this case saxions and axinos are 
mainly thermally produced and are short-lived, thus suppressing their contributions 
to AN e ff and fi^/i 2 . Once f a > 10 12 GeV, axinos become long-lived and enhance the 
neutralino abundance. In most cases, the augmentation leads to an overabundance of 
neutralinos. However, for very high values of f a (> 10 15 GeV) and small values of sq 
(< 10 -3 ), both the thermal production of saxions and axinos and the production of 
saxions via coherent oscillations become suppressed, thus resulting effectively in the 
usual thermal production of neutralinos accompanied by CO-produced axions. Then, 
once again, the DM, BBN and AN e ff constraints can be simultaneously satisfied. If 
&z h 2 is less than the observed dark matter relic abundance, the remaining DM 
abundance may be comprised of CO-produced axions, once the appropriate value of 
the misalignment angle 6i is chosen. These SUA SUSY models tend to have large 
enough direct /indirect WIMP detection rates that they should be soon seen by such 
experiments. 

Finally, we point out that if an excess on N e tf is confirmed by Planck data, it can 
be explained by a significant production of saxions via coherent oscillations followed by 
s — > aa decays. However, saxion decays to axinos must be suppressed in order to avoid 
subsequent entropy injection from axino decays. This is naturally satisfied if m s < 2m a . 

Note added: After completion of this work, we noticed Ref. |5^] appeared on a similar 
topic. 



- 25 - 



Acknowledgments 

This research was supported in part by the U.S. Department of Energy and FAPESP. 



A. Boltzmann equations for number and energy densities 

The general Boltzmann equation for the number distribution of a particle species can be 
written as j54j (assuming isotropy): 

B F BF 

^-H P ^ = Cl [F t ,F 3 ,p] (A.l) 

where Fi(p) is the number distribution of particle i as function of momentum p, C represents 
a source/sink term and H is the Hubble constant: 



1 pr_ 



H=J-^ T (A.2) 
V 3 Ml 



with px = YliPi- The number, energy and pressure densities are given in terms of as: 

dp 
2^ 



Pi{t) = j ^EiFifat) (A.3) 



where m« is the mass of particle i and Ei = ypf + mf . Using Eq. ( |A.1| ) we obtain the 
following equations for the number and energy densities: 

ir +3Hn * = J 2^ pCl 

d Pi , orr/_ , r>\ f d P 2 



dt + 3H( Pi + P l ) = J ^p'EiQ. (A.4) 
We will assume that C is given by: 

C = Cdec + C co ll (A-5) 

where Cdec contains the contributions from decays (i — > j + X and j — > i + X) and C co u 

from collisions with the thermal plasma. Below we compute each term separately, under 
some simplifying assumptions. 

A.l Collision Term 



The collision term C co ii for the i + j -<->• a + b process is given by |54| : 

/ ^P 2<J coii = - J dUidU.dU.dn^TT^S^ip, + Pj - Pa - Pb )\M\ 2 

x [FiFjil ± F a )(l ± F b ) - F a F b (l ± F t )(l ± Fj)] (A.6) 
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where dU { = d 3 pi/((2vr) 3 2 J E i ). Since we are ultimately interested in Eqs. (|A.4j) for the 
number and energy densities, we will consider the following integral: 

/ Z^ p2Cco11 X E ? =-I du i du J du adU b (2TT) i 5W( Pi +p j - Pa - Pb )\M\ 2 

x [FiFj(l ± F a )(l ± F b ) - F a F b (l ± ± Fj)] x £f (A.7) 

where a = 0(1) for the contribution to the number (energy) density equation and the plus 
(minus) sign is for bosons (fermions). Below we derive a simplified version of the above 
equation, valid under some approximations. The first approximation assumes F{ <C 1, so 
1 =t Fi ~ 1. Furthermore, here we assume that the distributions can be approximated by : 

Fi~exp(-(Ei-Hi)/T) (A.8) 

so the collision term can then be written as: 

/ ^2P 2 CcoiiET = - (exp((^ + N )/T) - exp((/i a + /i 6 )/T)) 



x J dUidnjdnadU^n)^^ ^ + Pj - Pa - p b )\M\ 2 exp(-(Ei + E 3 )/T) x Ef 

where above we have used conservation of energy (E^ + Ej = E a + E b ). Since for the cases 
of interest the equilibrium distributions have zero chemical potential, we have: 



71 ■ 

^ = expim/T) (A.9) 
m 



so: 



2ir z \riinj n a n b J 

x J dn i dO j dn a dn h (2ir) 4 S^ ( Pi + Pj -p a - Pb )\M\ 2 exp(-(Ei + Ej)/T) x E, 



In particular, for the process i + i a + b, where a and b are in thermal equilibrium 
(jl a = Hb = 0): 

/^o-^ = -(|-i) 

x J dliid^dliadli^fS^ {pi + Pj -p a - Pb )\M\ 2 exp(-(Ei + E 3 )/T) x Ef 

= -(n 2 -n 2 )(cvE«). (A.10) 

For a = 0, the above equation is the well known contribution from thermal scatterings to 
the collision term. To estimate its value for a = 1, we assume: 

(avE) ~ (<7«)(£i> = (<tv)^ (A.ll) 
where ( ) represents thermal average. Thus: 

Cco ^ " ^ " lav)* , for a = 1 • (A " 12) 



7 This approximation is only valid for particles with a thermal distribution. However, since the collision 
term is responsible for keeping the particle i in thermal equilibrium with the plasma, it is reasonable to 
assume a thermal distribution for i while the collision term is relevant. 
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A. 2 Decay Term 



Now we derive a simplified expression for the decay term, under approximations similar to 
the ones used in the last section. The decay term includes the contributions from particle 



decay and injection from other decays and is given by [30]: 



( ,( />■ / ) - y/^I-'i + BR(j i)Tj / dqFj ( x /{p + q) 2 - mj, t 



rnj 

Jm 2 /4p 



(A.13) 



where BR(j — > i) is the branching ratio for j — > i + X decay times the multiplicity of i 
particles in the final state 8 . Once again we consider the integral: 



+BR(j -> i)T jmj f ^Ef J 2 ^ dqFj [^(p + q. 



2 — m 2 , t 



(A.14) 



with a = 0(1) for the contribution to the number (energy) density equation. The first term 
in Eq. (A.14) gives: 



2^ P R MP,t) ' 



-TiUiiriii^-) , for a = 
—TimiTii , for a = 1 



(A.15) 



The second term in Eq. ( A.14 ) can be simplified if we note that, for rrij 3> mj, the i 
particles injected from j decays are relativistic, so we have Ef ~ p a . The integrals over p 
and q can be conveniently rewritten through the following change of variables: 



= JP 2 + 



p and p = 



p2 + 



rn 7 



Pcos I 



After performing the integration over cos#, we obtain: 



(A.16) 



—p a 

2-7T 2 !,,.• 



rn 2 /4p 



dqFj I J (p + q) 2 - mj, t 



nj(-g-) , for a = 
rij/2 , for a = 1 



(A.17) 



Finally, replacing Eqs. ( A.15| ) and (|A.17 ) in Eq. ( A.14 ) and assuming (l/E) ~ 1/(E) = n/p 
(as in Eq. ( A.ll|) ), we have: 



^p 2 C dec {p,t)Ef 



-Tirmnf/pi + BR(j -)■ i)T jnijU 2 / pj , for a = ^ 
—TimiUi + BR(j — > i)Y jirnjUj /2 , for a = 1 



8 Eq. (A.13) assumes a two body decay kinematics of the type j — > j + X , where X can be equal to i 
but m,i,rax <C frij. Decays with small mass splitings or three body decays are not described by Eq. (A.13) 
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A. 3 Number and Energy Density Equations 

Using the results of Eqs. ( |A,12| ) and ( |A.18| ) in the Boltzmann equations for m and pi 
(Eq. flO))), we obtain: 

+ 3Hm = (nf - nfj (<rv) - T l m i ^- + V BR(j -> i)T jmj ^ (A. 19) 

Pi . , . Pj 

3+t 



*£L + 3 # (p . + P .) = _ ^ M Pi _ r . m . n . + ^ - )r 

3+1 



It is convenient to use the above results to obtain a simpler equation for pi/rii'. 

Besides the above equations, it is useful to consider the evolution equation for entropy: 

2tt 2 

S^—g* s (T)T 3 R 3 (A.21) 
45 

where R is the scale factor. With the above definition we have [54|: 

s= (^ 5 (r)|y /3 i? 4 ^i?(^x)lr^ 

^S = — Y,R{i^ x ) v i^i (A.22) 

i 

where R(i —> X) is the fraction of energy injected in the thermal bath from i decays. 
Defining: 

x = ln(R/R ), Ni = ln(nj/s ), and N s = ln(5/5 ) (A.23) 
we can write Eqs. ( |A~22l) , (|A~T9|) and (|A~20| ) as: 

N 's = J[fY. R{i ^ X ) T * m i ex P[^ + 3 * " ^s] ( A - 24 ) 

^ = -3- + Ef«(i - + *Y - 1] (A.26) 

H Pi/ni *—f. Hpj/njrii H \ ni J 

3t z ' 1 ' 

tf. = -3^ + V -»■ ij^m^ f J - (A.26) 

3+i J / 

where ' = d/dx. 

The above equation for N{ also applies for coherent oscillating fields, if we define: 

Ni = ln(nj/s ), and n { = pi/rm (A. 27) 

so 

R'i = (A. 28) 
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where we assume that the coherent oscillating component does not couple to any of the 
other fields. 

Collecting Eqs. flO^ - fljOtj) and flOg ), we have a closed set of first order differential 
equations: 



Entropy: 



N' s = — R{* -> X)T imi exp[Ni + 3x- N s ] 



(A.29) 



Thermal fields: 
Nj = -3 



r, m; , v^nn/- -\ r j m i n i , ( av )i 



H pi/n 



H Pj/rij rii H 



R' 



' H ... 

3TI 



Coherent Oscillating fields: 



■\ T j n i 
H m 



nj_Pi_ 

Pj Tli 



N' -- 
R' = 0. 



3 -i 



5| -il 



(A.30) 



(A.31) 



As seen above, the equation for Ri = pi/rii depends on Pi/m. A proper evaluation of 
this quantity requires knowledge of the distribution Fi(p,t). However, for relativistic (or 
massless) particles we have Pi = pi/3, as seen from Eq. ( |A.3j ), whilst for particles at rest 
we have Pi = 0. Hence Fi(p,t) is only required to evaluate the relativistic/non-relativistic 
transition, which corresponds to a relatively small part of the evolution history of particle 
i. Nonetheless, to model this transition we approximate Fi by a thermal distribution and 
take Ti, pi <C rrii, where is the temperature of the particle (which can be different from 
the thermal bath's). Under these approximations we have: 



— = Tj and — 

m rii 



K^rrii/Ti) rrii 



+ 3 



K 2 {mi/Ti) Ti 

where K\i are the modified Bessel functions. In particular, if rrii/Ti 3> 1: 



EL 

m 



Ti 



3 rrii 
- + ^ + 3 



P 



2rrii ( R 



I 



(A.32) 



(A.33) 



rii 3 \ in; 

As shown above, for a given value of Ri = Pi/rii, Eq. QA.32 ) can be inverted to compute 
Ti (= Pi/m): 

P (A.34) 



Pi 

— = T{Ri 
rii 



Since we are interested in the non-relativistic/relativistic transition, we can expand the 
above expression around Ri/rrii = 1, so Pi /rii can be written as: 



rii 



2rrii ( R 



— - 1 ) + rrii V] a n ( — 
rrii I ' \rrii 



71>1 



Ri 



(A.35) 
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where the coefficients a n can be numerically computed from Eq. ( |A,32 ). The above approx- 
imation should be valid for rrii/Ti > 1 (or Ri > wij). On the other hand, for mj/Tj <C 1 
(or Ri 3> mj), we have the relativistic regime, with Pi/ Hi = Ri/3. Therefore we can 
approximate the Pi/rii function for all values of Ri by: 



£ = J ¥(S- 1 )+^En>l«n(g-l) j0T Rl <R 

n i \ , for Ri > R 



(A.36) 



where the coefficients a n are given by the numerical fit of Eq. ( A.32| ) and R is given by the 
matching of the two solutions. 

Finally, to solve Eqs. ( |A.29| )-( |A.31| ) we need to compute H according to Eq. ( |A,2j ), 
which requires knowledge of the energy densities for all particles (pi) and for the thermal 
bath (pr). The former are directly obtained from Ni and Ri, while the latter can be 
computed from Ng: 

// " ,( / /,>) * 3 T R eMNs/3 -x]^ PR = ^ff*(T)T 4 . (A.37) 



g*s(T)J 1 ,IX 30* 

Eqs. (|A.29|) - (|A.31D , with the auxiliary equations for H (Eq. flAjp ) and Pi/rii (Eq. (|A.36|) ) 
form a set of closed equations, which can be solved once the initial conditions for the num- 
ber density (nj), energy density (pi) and entropy (S) are given. For thermal fluids we 
assume: 

m(T R ) = {\ ' if M«*/2^ < 2 {A38) 

{ ni(T R ) , if (av)ini/H\ T =T R > 2 
— (Tr) = ^(T R ) (A.39) 

where is the equilibrium energy density (with zero chemical potential) for the particle i. 
For coherent oscillating fluids, the initial condition is set at the beginning of oscillations: 

nW) = ^ (A.40) 
^(T t osc ) = ra, (A.41) 

m 

where T° sc is the oscillation temperature, given by 3H(T° SC ) = rrii(T° sc ) and pf the initial 
energy density for oscillations. 

Finally, the initial condition for the entropy S is trivially obtained, once we assume a 
radiation dominated universe at T = Tr: 

2tt 2 

S(Tr) = —g,(T R )T*B*. (A.42) 
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